Tarea 2

Tarea 2: Frequentist Inference

CC6104: Statistical Thinking

Integrantes :

• Franco Cattani
• Juan Molina

Cuerpo Docente:

• Profesor: Felipe Bravo M.
• Auxiliares: Martín Paredes, Benjamín Farías
• Ayudantes: Renzo Zanca, Kevin Iturra, Emilio Díaz, Scarleth Betancurt

Fecha límite de entrega:

Índice:

1. Objetivo
2. Instrucciones
3. Referencias
4. Elaboración de Código

Objetivo

Bienvenid@s a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

• La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacerla de forma individual.
• La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
• El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
• Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
• No serán revisadas tareas desarrolladas en Python.
• Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
• Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.

Referencias:

Slides de las clases:

• Introduction to Statistical Inference
• Design of Experiments & Hypothesis Testing

Enlaces a videos de las clases:

• Introduction to Statistical Inference: video1 video2 video3 video4 video5 video6 video7 video8
• Design of Experiments & Hypothesis Testing: video1 video2 video3 video4 video5 video6

Documentación:

• Boostrap en r
• ggplot2

Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:

# Manipulación de estructuras
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# Para realizar plots
library(ggplot2)
library(plotly)

## 
## Adjuntando el paquete: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

# Manipulación de varios plots en una imagen.
library(gridExtra)

## 
## Adjuntando el paquete: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

Elimine eval=FALSE para ejecutar las celdas.

Pregunta 1: Estimadores.

• Considere una serie de variables aleatorias \(X_i\) que siguen una distribución de Bernoulli de parámetro \(p=0.5\) y un estimador de \(p\) dado por \(\hat{p}_{n, \sigma} = \epsilon_{\sigma} + \frac{1}{n} \displaystyle{\sum_{i=1}^{n}}X_{i}\) donde \(\epsilon_{\sigma} \sim \mathcal{N}(0,\sigma)\). Para \(\sigma = 0, 1, 2, 4\) grafique el valor de \(\hat{p}_{n, \sigma}\) y comparelo con el valor verdadero. Calcule el sesgo del estimador para cada valor de \(\sigma\).

# Function to calculate bias
calculate_bias <- function(estimates, p_true = 0.5) {
  mean(estimates) - p_true
}

# sigma = 0
cte <- rep(0.5, 1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 0)
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_0 <- data.frame(x = x, p_hat = data)
bias_0 <- calculate_bias(data)

# Create individual plots
p1 <- ggplot(df_0, aes(x = x, y = p_hat)) +
  geom_line(color = "blue") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 0)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 1
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 1)  # sigma = 1
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_1 <- data.frame(x = x, p_hat = data)
bias_1 <- calculate_bias(data)

p2 <- ggplot(df_1, aes(x = x, y = p_hat)) +
  geom_line(color = "red") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 1)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 2
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 2)  # sigma = 2
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_2 <- data.frame(x = x, p_hat = data)
bias_2 <- calculate_bias(data)

p3 <- ggplot(df_2, aes(x = x, y = p_hat)) +
  geom_line(color = "green") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 2)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 4
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 4)  # sigma = 4
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_4 <- data.frame(x = x, p_hat = data)
bias_4 <- calculate_bias(data)

p4 <- ggplot(df_4, aes(x = x, y = p_hat)) +
  geom_line(color = "purple") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 4)", x = "n", y = "Estimator of p") +
  theme_minimal()

# Print plots
print(p1)

print(p2)

print(p3)

print(p4)

# Print out the biases for each sigma
cat("Bias for sigma = 0:", bias_0, "\n")

## Bias for sigma = 0: -0.0005238953

cat("Bias for sigma = 1:", bias_1, "\n")

## Bias for sigma = 1: -0.01096729

cat("Bias for sigma = 2:", bias_2, "\n")

## Bias for sigma = 2: 0.01166566

cat("Bias for sigma = 4:", bias_4, "\n")

## Bias for sigma = 4: 0.08243936

Respuesta

• ¿Qué puede decir del sesgo de los estimadores \(\hat{p}_{n,\sigma}\)? ¿Son los estimadores con menos sesgo más precisos? Justifique.

Respuesta

Estos muestran alta variabilidad, tanto para altos como bajos valores de \(\sigma\). El único que mantiene cierta regularidad es \(\sigma = 0\). Esto se debe a que para \(\sigma = 0\) no hay ruido en el estimador, por lo que es esperable que este se acerque al valor real a medida que aumenta n. Vemos que los otros valores de sigma si se diferencian del valor real, aumentando esta diferencia promedio a medida que aumenta n (para \(\sigma = 1\) y \(\sigma = 2\) es cercana a \(\pm 4\), mientras que para \(\sigma = 4\) llega hasta \(\pm 15\)). En conclusión, los estimadores con menos sesgo si son más precisos.

Pregunta 2: Intervalo de Confidencia

En las preguntas 2, 3 y 4 deberá trabajar con el dataset diabetes_prediction_dataset.

datos <- read.csv("diabetes_prediction_dataset.csv", header = TRUE)
head(datos)

Una forma de estimar la distribución de la media de una población es a través de la realización de la sampling distribution de una distribución cualquiera. El valor obtenido en cada sampling distribution nos entrega un estimador de la media que posee una determinada distribución de la población. Sabiendo esto, calcule la sampling distribution de las columnas bmi y blood_glucose_level, obteniendo el intervalo de confianza de \(95\%\) para cada una de las medias obtenidas desde la distribución. Para realizar este experimento considere los siguientes puntos:

• Obtener las medias y desviaciones estándar de la población real.
• Realizar una sampling distribution con un tamaño de muestra igual a \(100\). Repita la obtención de la media un número elevado de veces (recomendación \(5000\) veces).
• Calcular el intervalo de confianza para cada uno de las medias obtenidas en las iteraciones.
• De acuerdo a los valores obtenidos (medias e intervalos de confianza), grafique una muestra aleatoria de las medias obtenidas en conjunto a sus intervalos de confianza. Aquí debe notar que, si el intervalo de confianza contiene a la media de la población, este se considerará como parte del intervalo de confianza del \(95\%\), haga un conteo de cuantos valores están contenidos en este intervalo. Para graficar los intervalos deconfianza puede ser útil utilizar forest plot.
• Compare las estimaciones para bmi y blood_glucose_level y señale que diferencias y similitudes encuentra entre estas (en cuanto a la calidad de la estimación, no los valores obtenidos). De una explicación de por qué cree que se dan las similitudes/diferencias.

Hints:

• Para comparar las 2 estimaciones puede ser útil normalizar bmi y blood_glucose_level a una misma escala.
• Para realizar la sampling distribution podría serle útil el comando sample().
• Puede ser útil generar un dataframe para graficar con ggplot2.

Respuesta:

# Definimos tamaños de muestreo
tamano_muestra = 100
n_muestras = 5000

# Inicializamos estructuras necesarias
list_mean_bmi = vector()
list_interval_bmi_left = vector()
list_interval_bmi_right = vector()
sucess_bmi = 0

list_mean_bgl = vector()
list_interval_bgl_left = vector()
list_interval_bgl_right = vector()
sucess_bgl = 0

# Obtenemos la media y desviación estándar de la población real
real_mean_bmi <- mean(datos$bmi)
real_sd_bmi <- sd(datos$bmi)

real_mean_bgl <- mean(datos$blood_glucose_level)
real_sd_bgl <- sd(datos$blood_glucose_level)

# Sampling distribution, calculo del intervalo de confianza y proporción.
for (i in 1:n_muestras) {
  
  # Generating samples
  sample_bmi <- sample(datos$bmi, tamano_muestra)
  sample_bgl <- sample(datos$blood_glucose_level, tamano_muestra)
  
  # Sampling distribution of the mean
  list_mean_bmi <- c(list_mean_bmi, mean(sample_bmi))
  list_mean_bgl <- c(list_mean_bgl, mean(sample_bgl))
  
  # Standard deviations of the samples
  sd_bmi <- sd(sample_bmi)
  sd_bgl <- sd(sample_bgl)
  
  # Calculating bmi confidence interval
  margin_error_bmi <- 1.96 * (sd_bmi / sqrt(tamano_muestra))
  ci_lower_bmi <- mean(sample_bmi) - margin_error_bmi
  ci_upper_bmi <- mean(sample_bmi) + margin_error_bmi
  list_interval_bmi_left <- rbind(list_interval_bmi_left, ci_lower_bmi)
  list_interval_bmi_right <- rbind(list_interval_bmi_right, ci_upper_bmi) 
  
  # Calculating bgl confidence interval
  margin_error_bgl <- 1.96 * (sd_bgl / sqrt(tamano_muestra))
  ci_lower_bgl <- mean(sample_bgl) - margin_error_bgl
  ci_upper_bgl <- mean(sample_bgl) + margin_error_bgl
  list_interval_bgl_left <- rbind(list_interval_bgl_left, ci_lower_bgl)
  list_interval_bgl_right <- rbind(list_interval_bgl_right, ci_upper_bgl)
  
  # Verify if the real mean is within the confidence interval
  if (ci_lower_bmi <= real_mean_bmi && ci_upper_bmi >= real_mean_bmi) {
    sucess_bmi <- sucess_bmi + 1
  }
  
  if (ci_lower_bgl <= real_mean_bgl && ci_upper_bgl >= real_mean_bgl) {
    sucess_bgl <- sucess_bgl + 1
  }
}

# Proporción de intervalos que contienen la media real
prop_bmi <- sucess_bmi / n_muestras
prop_bgl <- sucess_bgl / n_muestras

cat("Proporción de intervalos que contienen la media real (BMI):", prop_bmi, "\n")

## Proporción de intervalos que contienen la media real (BMI): 0.9474

cat("Proporción de intervalos que contienen la media real (Blood Glucose Level):", prop_bgl, "\n")

## Proporción de intervalos que contienen la media real (Blood Glucose Level): 0.953

# Convertir los datos a dataframes
df_bmi <- data.frame(
  mean = list_mean_bmi,
  lower = list_interval_bmi_left,
  upper = list_interval_bmi_right
)

df_bgl <- data.frame(
  mean = list_mean_bgl,
  lower = list_interval_bgl_left,
  upper = list_interval_bgl_right
)

df_bmi <- head(df_bmi, 100)
df_bgl <- head(df_bgl, 100)

# Columna que nos indica si la media real está dentro del intervalo de confianza calculado
df_bmi$outside_ci <- ifelse(df_bmi$lower > real_mean_bmi | df_bmi$upper < real_mean_bmi, "Fuera", "Dentro")
df_bgl$outside_ci <- ifelse(df_bgl$lower > real_mean_bgl | df_bgl$upper < real_mean_bgl, "Fuera", "Dentro")

# Normalización de las columnas
# Esto no lo usé al final pero lo dejé igual
normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}

# Plot de Intervalos de confianza
# Gráfico para BMI
ggplot(df_bmi, aes(y = 1:100, x = mean)) +
  geom_point(aes(color = outside_ci)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper, color = outside_ci), height = 0.2) +
  scale_color_manual(values = c("Dentro" = "black", "Fuera" = "red")) +
  geom_vline(xintercept = real_mean_bmi, color = "red", linetype = "dashed") + 
  labs(title = "Intervalos de Confianza para BMI", x = "Valor", y = "Muestra") +
  theme_minimal()

# Plot de Intervalos de confianza
# Gráfico para Blood Glucose Level
ggplot(df_bgl, aes(y = 1:100, x = mean)) +
  geom_point(aes(color = outside_ci)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper, color = outside_ci), height = 0.2) +
  scale_color_manual(values = c("Dentro" = "black", "Fuera" = "red")) +
  geom_vline(xintercept = real_mean_bgl, color = "red", linetype = "dashed") +
  labs(title = "Intervalos de Confianza para Blood Glucose Level", x = "Valor", y = "Muestra") +
  theme_minimal()

Respuesta La proporción de intervalos que contienen la media real parece ser muy buena y similar para ambas variables (Cercana al 94%). En general, ambas variables parecen mostrar resultados proporcionalmente similares a pesar de estar en una escala disitinta. Ambas variables suelen mostrar una cantidad baja de intervalos que no contienen la media real, a pesar de estar en una escala de magnitud distinta. Esto se debe a que se están haciendo la misma cantidad de experimentos, y estos al ser muchos, tienden a mostrar una proporción positiva en cuanto a contener la media real dentro del intervalo de confianza.

Pregunta 3: Estimación de Máxima Verosimilitud

El objetivo de esta parte será obtener los parámetros de las distribuciones de la media y la mediana de bmi. Para responder la pregunta realice los siguientes puntos:

• Simule la distribución empírica a través de mustras con repetición, como las de la pregunta anterior, y grafique los histogramas de la media y mediana.
• Utilice la estimación logaritmica de máxima verosimilitud negativa para obtener los parámetros de estas distribuciones.
• Gráfique a traves de un gráfico de calor el rango de valores en que se mueve la solución del problema de likelihood (puede ser útil persp() o filled.contour()).
• Aplique el comando nlminb() sobre la likelihood y encuentre el máximo o mínimo del problema a optimizar.
• Genere muestras lo suficientemente grandes de las distribuciones normales con los parámetros obtenidos por máxima verosimilitud, y comparelas con la distribución empírica mediante histogramas.

Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo la librería MASS).

Respuesta

# Obtenemos métricas de muestreo con repetición
means <- c()
medis <- c()
n_muestras <- 10000
tamano_muestra <- 100

for (i in 1:n_muestras) {
  
}

# dataframe con las medias y medianas de las muestras

# Histograma media

# Histograma mediana

# Media
# función log likelihood
ll_plot <- function(mu, sigma) {
  # ...
}

ll_plot <- Vectorize(ll_plot)

# definir espacio donde se va a evaluar ll_plot
mu_m <- 
sigma_m <- 
likelihood_m <- outer(X=mu_m, Y=sigma_m, ll_plot)

filled.contour(x=mu_m, y=sigma_m, z=likelihood_m, xlab=expression(mu), 
               ylab=expression(sigma), color = topo.colors)

# Mediana
# función log likelihood
ll_plot <- function(mu, sigma) {
  # ...
}

ll_plot <- Vectorize(ll_plot)

# definir espacio donde se va a evaluar ll_plot
mu_d <-
sigma_d <-
likelihood_d <- outer(X=mu_d, Y=sigma_d, ll_plot)

filled.contour(x=mu_d, y=sigma_d, z=likelihood_d, xlab=expression(mu), 
               ylab=expression(sigma), color = topo.colors)

likelihood_mean <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # ...
}

# Agrrgue el rango donde operará nlminb
nlminb(objective=likelihood_mean, start=c(, ), lower=c(, ), upper=c(, ))

likelihood_med <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # ...
}

# Agrrgue el rango donde operará nlminb
nlminb(objective=likelihood_med, start=c(, ), lower=c(, ), upper=c(, ))

# Muestra de medias


# Histogrmas

# Muestra de medianas


# Histogrmas

Pregunta 4: Bootstrap I

En esta pregunta se error e inetrvalo de confianza para la varianza de la columna HbA1c_level.

Las actividades por realizar son las siguientes:

• Programar el método Bootstrap para calcular el error estándar de la varianza. No se permite la utilización de librerías de bootstrap para esta parte.
• Visualizar a través de un de puntos las distintas varianzas obtenidas al realizar el muestreo con Bootstrap y comparar con el valor real.
• Visualizar a través de un gráfico el histograma obtenido al realizar el muestreo con Bootstrap y comparar con el valor real.
• Obtener el 95% de intervalo de confianza de la estimación.
• ¿Qué puede se inferir a pertir de los resutados de Bootstrap?

Respuesta:

# Bootstrap
B <- 5000
largo <- length(datos$HbA1c_level)
output <- c()

for (it in 1:B) {
  
}

# Gráfico de puntos

# Histogrma

# Obtenemos error e intervalos de confianze
se_vars <-
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE de la varianzas es: (%.3f)', se_vars)

Respuesta

Pregunta 5: Bootstrap II

El siguiente código genera una regresión lineal de bmi con respecto a age, es decir, una función lineal de age, \(y(age) = b + m\cdot age\), que pretende determinar el valor de bmi.

linearMod <- lm(bmi ~ age, data=datos)
print(linearMod)

## 
## Call:
## lm(formula = bmi ~ age, data = datos)
## 
## Coefficients:
## (Intercept)          age  
##    23.15536      0.09945

m <- linearMod$coefficients["age"]
b <- linearMod$coefficients["(Intercept)"]

ggplot() +
  geom_point(aes(x=datos$age, y=datos$bmi)) +
  geom_line(aes(x=datos$age, y=datos$age*m+b), color="red") +
  ggtitle("Regresión lineal") +
  ylab("bmi") + 
  xlab("Age")

Repita el proceso de la pregunta 4 para estimar el error e intervalos de confianza de los parámetros de la regresión (\(m\) y \(b\)). ¿Qué indican los resultados de Bootstrap?¿Un valor bajo en el error significa que la regresión es buena?

# Bootstrap
B <- 500
largo <- 1000
output_m <- c()
output_b <- c()

for (it in 1:B) {
  
}

# Obtenemos error e intervalos de confianze
se_m <- 
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE de la varianzas es: (%.3f)', se_m)

# Obtenemos error e intervalos de confianze
se_b <- 
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de los sesgos es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE del sesgo es: (%.3f)', se_b)

Respuesta

 

---

A work by CC6104